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1 .  Introduction. 

Two  point  boundary  value  problems  are  often  solved  by  finite 
difference  or  finite  element  methods.  Sometimes  one  can  find  a 
"closure"  or  "limit"  of  the  method  as  the  mesh  size  h — >0.  Such 
closures  are  often  shooting  methods  which  provide  more  flexibility 
and  adaptability  than  the  finite  difference  or  finite  element 
method.  The  double  sweep  method  (see,  for  example,  Babuska, 

Prager  and  Vitasek  [3],  Babuska  amd  Majer  [21,  Keller  and  Lentini 
[10] )  is  an  example  of  this,  which  is  also  related  to  the  invar¬ 
iant  imbedding  method  (see  Scott  [11]). 

In  this  paper,  we  shall  give  a  closure  of  the  Sturm  sequence 

V*t 

algorithm  for  computing  the  ntn  eigenvalue  and  eigenfunction  of 
the  finite  difference  discretization  of  a  Sturm-Liouville  problem. 
The  resulting  shooting  method  is  related  to  the  critical  (or  char¬ 
acteristic)  lengths  in  the  invariant  imbedding  method.  However,  a 
simple  counting  of  critical  lengths  does  not  produce  a  correct 
algorithm  for  calculating  the  nth  eigenvalue  (for  general  boun¬ 
dary  conditions).  We  shall  discuss  this  in  detail  in  section  3. 

In  addition,  our  method  gives  an  a-posteriori  error  estimate  for 
the  approximate  value  of  the  n*'1  eigenvalue. 

Standard  methods  by  finite  differences  or  finite  elements 
usually  find  the  first  n  eigenvalues  tj.  ,  1  ^  k  n.  Shooting 
methods,  such  as  those  discussed  in  Keller  [9],  find  one  eigen¬ 
value,  but  do  not  determine  which  one  it  is.  We  are  aware  of  only 
one  other  shooting  method  which  finds  the  nth  eigenvalue.  That 
method,  which  uses  the  Prufer  transformation,  was  given  by  Bailey 
[5]  and  Godart  [7],  and  was  also  discussed  by  Scott,  Shampine  and 
Wing  [12].  Our  method  has  the  advantage  that  it  generalizes  to 
higher  order  problems  and  systems.  Such  generalizations  will  be 
discussed  elsewhere. 

In  section  2  of  this  paper,  we  formulate  the  problem  and 
describe  the  shooting  method  (which  is  a  closure  of  the  Sturm 


sequence  algorithm).  In  section  3,  a  plausible  but  incorrect 
algorithm  is  discussed.  This  is  related  to  critical  lengths.  In 
section  4,  we  discuss  closures  of  numerical  processes.  In  this 
context,  we  outline  the  proof  of  our  main  theorem  (Theorem  1), 
which  is  the  basis  for  our  shooting  method.  In  sections  5  and  6, 
we  recall  details  of  the  finite  difference  discretization  and  the 
Sturm  sequence  algorithm.  Section  7  deals  with  the  closure  of  the 
Sturm  sequence  algorithm,  and  the  proof  of  Theorem  1.  In  section 
8,  we  summarize  our  results,  and  state  some  conclusions. 
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A  Shooting  Method  for  Eigenvalues  and  Eigenfunctions 
Consider  the  eigenvalue  problem 


(2,1a)  -(p(x)u')'  +  q(x)u  =  \  r(x)u  ,  for  0  <  x  1. 

(2.1b) 


aQ  u(0)  +  aQ  u'  (0)  =  0, 


aj  u ( 1 )  +  u' ( 1 )  =  0 , 
where  we  suppose  that: 

(a)  p(x) ,  q(x)  ,  r(x)  are  continuous,  and  p(x)  is  piecewise 
differentiable  on  [0,1], 

(b)  p(x)  >  0  and  r(x)  >  0  on  [0,1]. 

(c)  The  boundary  conditions  are  normalized  as  follows:  If 
3^  *  0,  then  3 1  =  1;  if  :i ^  -  0 ,  then  =  1 

( for  i  =  0,1). 

It  is  well-known  that  the  eigenvalues  of  (2.1)  form  an  increasing 
sequence : 


X1  <  *2  <  V3  < 


<  xn  < 


and  lim  v  , 

n->«  n 


Remark  1 ■  Without  any  difficulty,  we  could  assume  that  p,  q,  r 
are  piecewise  continuous  functions,  instead  of  continuous  func¬ 
tions.  For  simplicity,  we  have  assumed  they  are  continuous. 


Remark  2 ■  The  boundary  conditions  in  (2. 
in  the  form 


(2.1c) 


aQ  u(0)  +  3q  p ( 0 ) u' ( 0 ) 
a1  u( 1 )  +  P( 1 )u' ( 1 ) 


could  also  be  written 


0, 

0, 


which  would  be  more  natural,  for  physical  reasons.  Obviously, 
conditions  (2.1b)  and  (2.1c)  are  equivalent.  We  have  chosen 
(2.1b)  for  simplicity. 

We  shall  now  describe  the  shooting  method  for  finding  the 
nth  eigenvalue  and  eigenfunction  of  (2.1).  Choose  a  number  Vq, 
and  let  un(x)  be  a  nontrivial  solution  of  the  initial  value 


problem : 


-(p(x)u' )'  +  q(x)u  =  Xr  r (x)u  ,  for  0  .  x  <  1 , 

(22)  ^ 

[aQ  u(0)  +  ?Q  u' (0)  =0. 

{Choose  any  Initial  conditions  which  satisfy  the  boundary  condi¬ 
tion  aQ  u(0)  +  :3q  u'(0)  =  0.)  Let  NQ  =  Nq(Aq)  be  the  number  of 

zeros  of  uQ(x)  in  the  open  interval  (0,1).  (Note  that  the 

uniqueness  theorem  implies  that  all  zeros  have  multiplicity  1 . ) 

Let  v  =  v ( A Q )  =ax  u0(l)  +  B1  u^(l),  and 

fo  if  vun{l)>0  or  v  =  0 , 
a  =  o(XQ)  =  {  0 

Jl  if  v  Uq(1)  .  0  and  v  *  0. 

(Note  that  Aq  is  an  eigenvalue  if  and  only  if  v  =  0.)  Finally, 

let  N(\q)  =  N0  +  a .  Thus,  N(\Q)  is  the  number  of  zeros  of 

u0(x)  in  (0,1),  with  a  correction  which  depends  on  the  boundary 
condition  at  x  =  1.  The  shooting  method  depends  on  the  following 
theorem,  which  will  be  proved  in  §7. 

Theorem  1 .  N(\q)  equals  the  number  of  eigenvalues  of  (2.1)  which 

are  less  than  \Q. 

We  can  now  describe  the  shooting  method: 

Step  0.  Find  values  LQ  <  RQ ,  such  that  N ( Lq )  =  n-1  and 
N ( Rq )  =  n.  This  implies  that  LQ  :  '  <  RQ. 

Step  k.  For  given  values  <  Rk_1(  with  N(Lk_j)  =  n-1  , 

N  ( Rj^-i )  =  n,  find  values  Lk,  Rj_,  such  that  N(Lk)  * 
n-1  ,  N(Rjc)  -  n  ,  L^^  .  L^  <  R^  :  and 

Rk-Lk  <  Rk-l“Lk-l • 

STOP  when  R^-L^  .  r,  where  r  is  a  given  tolerance. 

To  implement  the  above  steps,  we  need  an  initial  value  solver 
to  compute  Uq(x)  and  a  nonlinear  solver  for  Lj.  and  R^.  The 
approximate  eigenvalue  An  will  be  either  the  midpoint  of  the 
last  interval  [L^.Ru],  or  the  last  approximation  v  found  by  a 


nonlinear  solver.  The  approximate  eigenfunction  <in(x)  will  be 
the  solution  of  (2.2),  with  ,\0  =  An. 

STEP  0  can  be  carried  out  either  by  using  estimates  of  *n 
which  the  user  may  have,  or  by  using  a  related  boundary  value 
problem  whose  eigenvalues  are  known.  For  STEP  k,  we  can  use  the 
bisection  method  for  the  integer-valued  function  N(\),  or  we  can 
combine  this  with  a  nonlinear  solver,  such  as  the  secant  method, 
for  the  function  v(\). 

Assuming  that  the  ODE's  under  consideration  are  solved 
exactly,  the  method  gives  a  sharp  a-posteriori  error  esinate  for 
the  eigenvalue  \n.  Of  course,  the  ODE's  will  be  solved  numeri¬ 
cally.  An  effective  implementation  of  this  method  must  relate  the 
accuracy  of  the  initial  value  solver  (governed  by  an  input 
tolerance  parameter)  to  the  value  R^-L^,  and  to  the  accuracy  of 
the  nonlinear  solver  for  finding  L^,  R^.  A  detailed  analysis  is 
required  to  determine  the  effects  of  the  approximate  solution  of 
ODE's  on  this  method,  and  on  the  a-posteriori  error  estimation. 
Such  an  analysis,  together  with  a  program  implementing  the  method, 
will  be  published  elsewhere. 

Remark  3 .  If  we  are  only  interested  in  the  eigenvalue  i  ,  and 
not  the  eigenfunction  on,  then  we  need  only  concern  ourselves 
with  the  count  of  the  zeros  of  u(x)  and  the  correction  term  a. 
This  can  be  obtained  by  solving  various  transformed  formulations 
(such  as  that  used  in  the  invariant  imbedding  method). 

Remark  4 ■  The  method  resembles  a  count  of  the  critical  lengths  in 
the  invariant  imbedding  method  (see  [11],  chap.  V).  However,  a 
simple  count  of  the  critical  lengths  (with  no  correction  term) 
does  not  produce  a  correct  algorithm,  for  general  boundary  condi¬ 
tions.  We  shall  return  to  this  point  in  's  3 . 

Remark  5 .  Theorem  1  leads  immediately  to  the  following  corollary, 
which  is  a  classical  theorem  about  eigenfunctions. 


J'JJUVVJ 


We  consider  an  alternative  algorithm  for  calculating  the  n^11 
eigenvalue  of  (2.1).  Let  uQ(x)  denote  a  solution  of  the  initial 
value  problem  (2.2),  as  before,  and  let  vQ(x)  =  u1uQ(x)  +  Uq(x). 
It  might  seem  natural  to  count  the  zeros  of  Vq(x),  rather  than 
counting  the  zeros  of  u0(x) ,  with  a  correction  using  v0(l),  as 
in  the  previous  section.  In  other  words,  if  N'(Vq)  denotes  the 
number  of  zeros  of  Vq(x)  in  (0,1),  we  might  conjecture  the 
following: 

Hypothetical  Theorem.  N'(\Q)  equals  the  number  of  eigenvalues  of 
(2.1)  which  are  less  than  \q. 

This  would  lead  to  an  alternative  algorithm  for  finding  *n,  by 
replacing  N(J  Q)  by  N'  ( A  0 )  STEP  0  and  STEP  k  (previous  sec¬ 

tion)  .  However,  the  above  Hypothetical  Theorem  turns  out  to  be 
false.  It  depends  on  monotonicity  properties  of  eigenvalues, 
which  are  valid  for  Dirichlet  boundary  conditions,  but  not  for 
general  boundary  conditions.  At  the  end  of  this  section,  we  shall 
give  an  example  where  the  alternative  algorithm  (using  N'Oq)) 
fails . 

Note  that  if  v0(x0)  =  0,  then  vQ  is  an  eigenvalue  on  the 
interval  [0,xQ],  and  xQ  is  a  critical  length  (as  used  in  the 
invariant  imbedding  method).  Thus  N'(\0)  As  the  number  of  cri¬ 
tical  lengths  (corresponding  to  A0)  which  are  less  than  1.  We 
see  that  a  simple  count  of  the  critical  lengths  does  not  produce  a 
correct  algorithm  for  calculating  the  nth  eigenvalue  of  (2.1). 

We  shall  now  indicate  the  role  played  by  monotonicity 
properties  of  eigenvalues.  For  0  <  y  1,  let  i  (y)  denote 
the  n  f  eigenvalue  of  the  same  operator  as  before,  with  the  same 
boundary  conditions,  but  on  the  interval  [0,y].  In  other  words, 

V  (y)  is  the  n**1  eigenvalue  for  the  problem: 


M 


-(p(x)u')'  +  q(x)u  =  A  r  (x)u  ,  for  0  <  x  j  y, 

(3.1)  aQ  u( 0 )  +  HQ  u' (0)  =  0, 
a1  u(y )  +  ,31  u'  (y)  =0. 

For  a  given  n,  we  ask  the  following  question: 

(Q)  Is  An(y)  a  decreasing  function  on  (0,1]? 

If  the  answer  to  (Q)  is  yes,  for  all  n,  then  the  Hypothetical 
Theorem  is  true  and  the  alternative  algorithm  is  correct.  For 
suppose  that  x0  (0  <  xQ  <  1)  is  a  zero  of 

vQ(x)  =  a j  u0(x)  +  :3 1  Ug(x),  where  uQ(x)  is  a  solution  of 

(2.2) .  Then  VQ  is  an  eigenvalue  on  [0,xQ],  i.e.,  \Q  =  Ak(xQ), 

for  some  k.  Since  A^ty)  is  a  decreasing  function, 

'k  =  <  lk^x0^  =  *0‘  In  t^is  waY>  each  zero  of  v0(x) 

corresponds  to  an  eigenvalue  less  than  \Q.  This  shows  that  the 
Hypothetical  Theorem  and  alternative  algorithm  are  correct,  if 
A n ( y )  is  a  decreasing  function,  for  all  n.  Unfortunately,  this 
is  not  true  for  general  boundary  conditions. 

If  the  boundary  condition  at  the  right  endpoint  is  a 
Dirichlet  condition:  u(y)  =  0,  then  the  classical  monotonicity 
theorem  tells  us  that  An(y)  is  a  decreasing  function,  for  all 
n.  But  Greenberg  [8]  has  shown  that  if  the  boundary  condition  at 
the  right  endpoint  is  not  a  Dirichlet  condition  (i.e.,  *  0), 

then  for  given  nQ  >  1,  there  exist  coefficient  functions  p(x) , 
q(x) ,  r(x)  and  a  subinterval  [a,b]  (0,2],  so  that  the  eigen¬ 
values  A.(y),  A9(y),  ...,  A  (y)  are  increasing  functions  in 

x  t  nQ 

[a,b] .  (On  the  other  hand,  for  given  p(x),  q(x) ,  r(x),  aQ,  vQ, 

Uj,  3 j ,  there  exists  n1  .  1,  so  that  for  n  ;  n1  ,  A n(y)  is  a 
decreasing  function  on  (0,1].)  Thus,  we  cannot  expect  the 
Hypothetical  Theorem  and  alternative  algorithm  to  be  correct  for 
general  boundary  conditions.  We  now  give  a  concrete  example  where 
they  fail. 
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Example . 

(3.2) 


For  0  <  y  <  1,  consider  the  eigenvalue  problem: 

-(p(x)u'  )'  =  X.  u  ,  for  0  <  x  y , 
u ( 0 )  +  u'  (0)  =  0  , 
u' (y)  =0. 


The  energy  norm  is  given  by 


(3.3) 

and 


B ( v , v) 


-p(0)v(0)2  + 


[  P(X) v'  ( x ) 2 dx / 
0 


(3.4) 


2 

where  v!:  = 

Putting  v(x) 


A.  :  ( y )  =  inf  , 

V^C1[0,y]  ,|V! 

y 

v(x) 2dx. 

0 

si,  we  find  that  A  ^  ( y )  - 


y 


Thus 


(3.5)  X 1 (y)  <  0. 

We  now  consider  the  two  algorithms  (given  by  Theorem  1  and  the 
Hypothetical  Theorem)  for  finding  the  number  of  eigenvalues  V  <  0 
(for  the  interval  [0,y]).  We  must  solve  the  initial  value 
problem : 

[-(p(x)u'  )  '  =  0 

(3.6) 

| u ( 0 )  =1  ,  u'  (0)  =  -1  . 


Denote  the  solution  by  u0(x),  and  let  vQ(x)  =  u'Q(x)  .  We 
obtain:  -p  u£,  =  constant  =  p(0),  so  that 


n 

X 

o 

> 

uo{x>  “  - 

(3.7) 

1 

X 

u0(x)  = 

1  -  1  pf?i  dt- 

0 

Since  Uq(x) 

< 

0  ,  u0(x)  is 

a  decreasing  function 

uQ  {  0  )  =  1 )  . 

For  a  given  y 

( 0  <  y  <  1 ) ,  either 

( A )  uQ ( y )  £ 

0 

and 

Ug (y )  < 

0 ,  or 

(B)  uQ(y)  < 

0 

and 

ufc(y)  < 

0 . 

1 1 


,»(  ,14  k|  «• 


>*  Jl  J»J  A*.W  I 


In  case  (A),  NQ  *  [number  of  zeros  of  uQ(x)  in  (0,y)]  =  0, 
a  =  1 ,  and  N  =  Nq  +  a  =  1 . 

In  case  (B),  NQ  =  1  ,  a  =  0,  and  N  =  NQ  +  a  -  1 . 

Thus  we  see  that  the  algorithm  of  Theorem  1  counts  1  negative 
eigenvalue  on  [0,y],  for  all  y  in  the  interval  (0,1]. 

On  the  other  hand,  the  alternative  algorithm,  based  on  the 
Hypothetical  Theorem  counts  the  zeros  of  vQ(x)  =  Ug(x)  in 
(0,y).  Since  vQ(x)  <  0  /  N'  =0,  predicting  no  negative  eigen¬ 
values!  Here  we  have  an  example  where  the  Hypothetical  Theorem 
and  alternative  algorithm  are  incorrect. 

Remark .  Greenberg  [8]  has  shown  that  in  the  above  example,  Ajfy) 
is  an  increasing  function  on  (0,1].  Thus  we  have  "reverse  mono¬ 
tonicity"  in  this  example! 
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4 .  Closure  of  a  Numerical  Process. 

Standard  methods  for  solving  boundary  value  problems  involve 
two  stages : 

(A)  A  discretization  method  (such  as  finite  differences  or  finite 
elements)  which  approximates  the  continuous  problem  by  a 
finite-dimensional  algebraic  problem. 

(B)  An  algorithm  for  solving  the  algebraic  problem. 

It  may  happen  that  the  algebraic  algorithm  is  equivalent  to 
finding  a  numerical  solution  u^  of  an  initial  value  problem,  and 
manipulating  this  solution  in  some  way  to  find  some  desired 
approximate  value  z^,.  Perhaps  z^  =  F^uh^'  where  is  some 

function.  As  the  mesh  size  h — >0,  suppose  that  Fh — >F  (some 
function),  u^ — »u  (the  solution  of  the  initial  value  problem), 
and  z^ — >z  (the  value  we  wish  to  compute).  Of  course  we  then 
have  z  =  F(u).  Thus  we  obtain  an  alternate  method  for  calculat¬ 
ing  z,  which  we  call  a  closure  of  the  numerical  process.  This 
method  consists  in  solving  the  initial  value  problem  to  find  u, 
and  then  calculating  F(u)  to  find  z. 

We  have  in  mind  situations  where  the  initial  value  problem  is 
hidden  in  the  numerical  process.  Finding  a  closure  amounts  to 
discovering  a  hidden  initial  value  problem.  If  we  find  a  closure, 
then  we  gain  a  better  understanding  of  various  properties  of  the 
algorithm,  and  we  can  improve  it.  The  algebraic  algorithm  may 
have  involved  a  low  order  method,  such  as  Euler's  method.  We  can 
choose  a  higher  order  initial  value  solver,  or  one  which  is 
especially  adapted  to  the  particular  problem. 

We  shall  show  that  Theorem  1  is  a  closure  of  a  numerical 
process.  In  the  next  section,  we  review  the  discretization  by 
finite  differences.  In  the  following  section,  we  discuss  the 
algebraic  algorithm,  using  Sturm  sequences.  In  this  context,  z^ 
is  the  number  of  eigenvalues  (of  the  finite  difference  equation) 
which  are  less  than  The  Sturm  sequence  turns  out  to  be 


related  to  the  initial  value  problem  (2.2).  counts  the  number 

of  sign  changes  in  the  Sturm  sequence.  This  turns  out  to  be  the 
same  as  the  number  of  sign  changes  in  the  numerical  solution  uh 
of  (2.2),  except  for  the  last  term.  This  last  term  gives  rise  to 
the  correction  term  a  (mentioned  in  §  2 ,  where  N(*0)  is 
defined) . 

The  general  notion  of  closure  is  probably  due  to  S.L.  Sobolev 
[13],  [14].  In  [14],  Sobolev  defines  closure  and  applies  it  to 
the  elimination  method  for  solving  a  discretized  integral 
equation.  Babuska,  Prager  and  Vitasek  [3],  [4]  find  closures  of 
several  algorithms.  The  closure  method  seems  to  have  many  other 
possibilities  for  application.  We  plan  to  return  to  some  of  them 
in  the  future. 


5 .  Discretization. 


We  shall  discretize  the  boundary  value  problem  (2.1)  by 
finite  differences.  We  use  uniform  mesh  of  size  h  =  1/n,  and 
the  finite  difference  operator: 

A  (V  =  £  (u  -  u  ,  ). 


1+7 


1_7 


Approximate  the  differential  equation  by  the  difference  equation 

-A  (pi  A  (ui)  +  q1ui  =  X  riui, 

or  equivalently: 

(5.1)  -p  +  (P.  L  +  P  L  +  h2qi)u1  -  p  Lui  +  1  =  \  h2riui. 

*  7  1~2  i+J  1+2 

The  boundary  condition  aQ  u(0)  +  3 Q  u'(0)  =0  is  approxi¬ 
mated  by  the  difference  equation 


a0  u0  +  '0 


ul"u0 


,?0U1 

=  0  '  or  u0  =  H0=K a~ 


(5.2) 


Substituting  this  into  (5.1),  with  i  =  1,  we  obtain 
-ha. 


3n-h an  Pt  +  P3  +  h  <*1  U1  '  P3  u2  “  x  h  rlul‘ 


The  boundary  condition  u(l)  +  B ^  u' (1)  =  0  is  approxi¬ 

mated  by  the  difference  equation 


«1  un  +  Bx 


'u  -u_  , 
n  n- 1  _ 

- K -  “ 


_  „  _ _  _  /Jlun-1 

-  0  ,  or  un  3  -fjjg  - 


(5.3) 


Substituting  this  into  (5.1),  with  i  =  n-l,  we  obtain 

ha 


1  2 

<P„  3  +  TTWT  i  +  h  ‘Jn-l,un-l 


_  s  n-2  _  ~  j  i  tiiw  -  ~  i 

n-j  n~7  1  1  n  7 


=  l  hzrn  ,u, 


n-l  n-l  ' 


Thus,  we  have  a  finite-dimensional  eigenvalue  problem: 
(5.4)  Au  =  A.  Ru , 

where  A  is  the  finite  difference  matrix 


‘J \0p  v«  !-<■ 


■ti-r 


with 

al  = 

so"htlo 

+  p3  +  h‘q1( 

2 

ai  " 

P,  i  +  P.  i  +  h2t?i  t 2 

H  1+T 

a 

=  n  4- 

ha, 

an-l 

P_  :»  + 

11  2 

'5l+hal  Pn-1  + 

bi  = 

-P<4  U 

5  i  <  n-2 )  , 

and 

R  is 

the  diagonal  matrix 

6.  Sturm  Sequences. 


The  finite  difference  discretization  has  approximated  the 
continuous  eigenvalue  problem  by  the  finite-dimensional  eigenvalue 
problem 

(6.1)  Au  =  A.  Ru. 

The  algebraic  algorithm  we  shall  use  for  this  is  the  method 
of  Sturm  sequences.  We  shall  briefly  recall  the  definition  and 
main  theorem  on  Sturm  sequences,  in  the  context  of  (6.1).  Details 
about  Sturm  sequences  can  be  found  in  Bathe,  Wilson  [6],  Stoer, 
Bulirsch  [15]  and  Wilkinson  [16]. 

For  a  given  number  AQ,  the  Sturm  sequence 

S0(AQ),  Sl(*0)'  Sn-l{X0)  is  deflned  bY:  S0(A0)  =  1,  and 

for  1  <  i  s,  n-1  ,  Sj(Aq)  is  the  leading  i  i  principal  minor  of 
A-.\q  R.  In  other  words,  S^(AQ)  *s  tbe  determinant : 


The  following  theorem  is  the  main  fact  that  we  shall  need. 

Sturm  Sequence  Theorem.  Let  c(vQ)  be  the  number  of  sign  changes 
in  the  sequence  S0(\Q),  S x ( A Q ) ,  ...,  sn_i(' q) •  Then  c ( a 0 ) 

equals  the  number  of  eigenvalues  (of  the  problem  (6.1))  which  are 
less  than 

Note:  If  si (\q)  =  0,  then  it  counts  as  a  sign  change,  in  the 

above  theorem. 


Closure  of  the  Sturm-Seauence  Alaorithm.  Proof  of  Theorem  1 


For  a  given  number  \0,  consider  the  initial  value  problem 


(7.1) 


-(p(x)u' )'  +  q(x)u  =  \q  r{x)u  ,  for  0 
u(0)  =  ,*0  ,  u'  (0)  =  -ct0 . 


Note  that  the  initial  conditions  were  chosen  to  satisfy  the 
boundary  condition  aQ  u(0)  +  ;3  Q  u'(0)  =  0.  Thus,  equation  (7.1) 
is  essentially  the  same  as  (2.2). 

Consider  a  finite  difference  approximation  to  (7.1).  We 
obtain  the  same  difference  equations  (5.1),  including  the  first 
equation  (5.2),  but  with  XQ  replacing  X.  We  do  not  obtain  the 
last  equation  (5.3),  because  u  is  not  required  to  satisfy  the 
boundary  condition  at  x  =  1.  Thus,  the  finite  difference 
approximation  satisfies: 


(7.2) 


(Uj-Xq  h^r, )un  +  b,u9  =  0, 


1 1 '  1 


(7.3)  bi-iui-i  +  (ai~*0  h2ri)ui  +  biui+i  =  0  ’  for  2  :  i  i  n-2 
It  also  satisfies  a  last  equation 

'Pn_|un-2  +  <Pn_3  +  Pn_L  +  h2<3n-l’un-l  '  *n_LUn 


=  h  r„  ,  u 


which  can  be  written 


f 

(  7  •  4  )  bn-2un-2  +  an-l  +  7^1^  Pn_L 


n-1  n-1 , 


-  *0  h  rn-l  un-l  +  bn-lun  = 


We  shall  relate  the  Sturm  sequence  Sq( < Q) , Sj ( 1 Q ),..., Sn_1 ( X 
to  the  solution  u^ ,  u2 ,  . . . ,  un  of  the  above  finite  difference 
equations.  We  first  consider  the  case  where  all  S^(Xq)  are 


nonzero.  We  shall  use  the  notation 


(7.5) 


Ai  =  ai  -  l0  h  rr 


For  1  •  i  x  n-2,  consider  the  first  i  equations: 
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’  s'  . 


then  (7.3)  shows  that  =  0.  Similarly,  this  would  imply 

that  u-  =  u2  =  . . .  =  =  0.  The  initial  conditions 

u(0)  =  3q,  u'  (0)  =  -cXq  in  (7.1)  are  implemented  in  the  approxi- 


u,-uQ 

mate  solution  by:  uQ  =  '?0  ,  ■  =  -«0-  Since  u1  =  0,  this 

implies  that  uQ  =  haQ  and  uQ  =  i3q,  so  that  r0  =  huQ.  This 
relation  can  be  satisfied  by  at  most  one  value  of  h,  since  uQ 
and  are  not  both  zero.  We  assume  that  h  is  small  enough  so 

that  J0  *  hciQ .  The  assumption  u^  =  ui  +  l  =  0  now  leads  to  a 
contradiction. 

We  have  shown  that  two  consecutive  terms  S^(Vq),  si+i('o^ 
cannot  both  be  zero,  nor  can  two  consecutive  terms  Uj,  ui  +  1> 

This  fact,  together  with  equation  (7.6),  written  as 


(7.14) 


uiSi<*0>  “  ~biui+lSi-l ( ' 0 : 


shows  that  S^.Vq)  =  0  if  and  only  if  ui  +  1  =  Now  suppose 

that  Sio(\Q)  =  Si^(VQ)  =  ...  =  Si  U0)  =0,  and  no  other  Si(VQ) 

is  zero.  The  argument  given  above,  for  the  case  where  all  S.jUq) 
are  nonzero,  shows  that  the  sequences 

sik+l'  Sik+2('o)'  Sik+1-1^0)'  and 

u,  Ui  . . . ,  u,  ,  have  the  same  number  of  sign  changes. 

1k+1  Ak+^  1k+l~1 


The  same  holds  for  the  beginning  and  end  segments  of  the  Sturm 
sequence,  corresponding  to  0  .  i  i  iQ-l  and  im  +  1  i  n-1. 
Thus,  if  we  count  a  zero  as  a  sign  change,  then  the  two  sequences 

s0(v0)'  sl('o>'  •••'  Sn-l(v0)  and  ul'  u2 . un-l'  vn  have 

the  same  number  of  sign  changes.  The  Sturm  Sequence  Theorem  now 
implies  that  the  number  of  sign  changes  in  the  sequence 
Uj ,  u2 ,  ....  un_j,  vr  equals  the  number  of  eigenvalues  of  the 

approximate  problem  Au  =  *  Ru,  which  are  less  than  aq. 

We  now  consider  what  happens  as  the  mesh  size  h — *0.  Let 
u0(x)  be  the  solution  of  (7.1),  and  v  =  u.  uQ ( 1 )  +  <1  Uq ( 1 ) .  We 
shall  consider  four  cases,  according  to  the  possibilities  that 


uQ(l)  and  v  are  zero  or  not. 


Case  1 .  Uq ( 1 }  *  0  ,  v  *  0 . 

The  continuous  problem  has  only  finitely  many  (say  m) 
eigenvalues  less  than  It  is  known  that  the  first  m  eigen¬ 

values  of  the  approximate  problem  converge  to  the  first  m  eigen 
values  of  the  continuous  problem  as  the  mesh  size  h — >0.  The 

number  of  sign  changes  in  the  sequence  u. ,  u^ . u__, 

converges  to  the  number  of  zeros  in  (0,1)  of  the  solution  u0(x 
of  the  initial  value  problem  (2.2)  and  (7.1). 

We  must  find  what  happens  to  vn  as  h — >0.  Recall  that 


vn 


13 1 

bn-lun  +  U 1  +ha  1  pn_i_un-l 


and  b 


calculation  shows  that 
(7.15) 


vn  = 


hpn-i- 


s  2+ha^ 


alun  +  11 1 


"-1  '  Pn-1-' 


u  -u 
n  n-1 


A  simple 


Recall  that  the  boundary  conditions  are  normalized:  either 
■3,-0  and  u1  =  l,  or  3 1  =  1  .  Thus 


vn  = 


n_r 


hpn-i 


l+h<Y. 


if 


=  0 


“lun 


’1 


=  1 


This  shows  that  for  small  h  ,  v  has  the  same  sign  as 
v  =  u0(l)  +  u^(l),  where  u g(x)  is  the  solution  of  (7.1). 

Thus,  for  small  h,  the  comparison  of  signs  between  un_i  and 
vn  is  the  same  as  that  between  uQ(l)  and  v.  This  shows  that 
the  correction  term  a  should  be  1  if  there  is  a  sign  change 
between  u(l)  and  v,  and  otherwise  t  should  be  0.  This 
concludes  the  proof  of  the  theorem,  in  Case  1. 


Case  2 .  uQ ( 1 )  -  0  ,  v  *  0. 

In  this  case,  o  =  1  by  definition, 
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1  and 


V  ”  ui  uo(1)  +  u0  ^ 1  ^  VX)  has  Posit:ive  values  in  sorr.e 

small  interval  (l->  ,1),  then  ( 1 )  <  0,  and 

v  =  (« 1  uQ  ( 1 )  +  Uq  ( 1 )  =  Uq  ( 1 )  <  0  .  If  uQ(x)  has  negative  values 
in  ( 1 —z  , 1 ) ,  then  Uq ( 1 )  >  0  and  v  >  0.  Thus  there  is  always  a 
change  of  sign  between  uQ(x),  for  x  close  to  1,  and  v. 

Now  consider  the  sequence  ,  u2 ,  un_1(  v  .  For  small 

h,  the  number  of  sign  changes  in  this  sequence  equals  the  number 
of  eigenvalues  (of  both  the  approximate  and  continuous  problems) 
which  are  less  than  \  q . 

(7.16)  #( eigenvalues  <  \q)  =  #(sign  changes  in  u1(u2 . un-l'vn^‘ 

Since  Uq ( 1 )  *  0  ,  Uq(x)  *  0  for  x  near  1.  Choose  x0 
near  1,  so  that  uQ(xc)  *  0  and  Uq(x)  *  0  for  xQ  x  i. 

Then  uQ(x)  is  monotone  on  [Xq,1]  ,  aiid  it  has  only  one  zero  in 
this  interval,  namely  at  x  =  1.  We  may  suppose  that  xQ  -  iQh. 
Then,  for  small  h. 


# (sign 

changes 

in 

ul- 

u2,  ... 

#( zeros 

of  uQ 

(x) 

in 

(0,x0)  ) 

#( zeros 

of  uQ 

(x) 

in 

(0,1)). 

Also,  from  (7.16),  we  have 


(7.18) 


#( eigenvalues  <  aq)  = 
#(sign  changes  in  u^ ,  u2 


#(sign  changes  in  u  =  ,  ..., 

1 0  “  1 


)  + 


Un_1-  Vn)> 


.,  un_,,  vn)  =  1. 


We  claim  that 

#(sign  changes  in  u,-  , 

0 

If  we  can  show  this,  then  (7.17)  and  (7.18)  imply  that  Theorem  1 
is  correct  in  this  case. 

Since  Uq(x)  *  0  on  [xQ,l],  the  sequence 


ui, 


+1 ,  ...,  un._j  is  monotone.  Thus  it  can  have  at  most  one 


sign  change.  If  there  is  a  sign  change,  then  either 
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(which  counts  as  a  sign  change)  or  un_i  has  the  opposite  sign  to 
that  of  Uq(x),  for  Xq  i  x  <  1  .  In  this  case,  there  is  no  sign 
change  between  un_^  and  vn  (since  we  have  seen  above  that 
there  is  always  a  sign  change  between  u(x),  for  x  near  1, 

1  ,  in 


and  v).  Thus  #(sign  changes  in  u;  , 

a0 


-i'  vn> 


'n-1 


this  case. 

If  there  is  no  sign  change  in  the  sequence 


u,  ,  u 


i^  +  1 


un-l'  then  un_1  and  v.  have  opposite  signs. 


Again  we  have  #(sign  changes  in  u(  .  u_  , ,  v_ )  =  1 .  We 

1  Q  il  i.  i  1 


have  shown  that 

#( eigenvalues  <  \0)  =  #(zeros  of  uQ(x)  in  (0,1))  +  1. 
This  proves  the  theorem  in  Case  2 . 


Case  3 .  uQ ( 1 )  *  0  ,  v  =  0 . 

In  this  case,  a  =  0  by  definition.  Also  is  an 

eigenvalue  of  the  continuous  problem,  say  \Q  =  >k.  Thus,  there 
are  k-1  eigenvalues  (of  the  continuous  problem)  less  than  vQ. 

We  must  show  that  the  eigenfunction  C>k(x)  =  uQ(x)  has  exactly 
k-1  zeros  in  (0,1).  Of  course,  this  is  a  classical  theorem. 

But  we  want  to  show  that  it  follows  from  Sturm  sequences. 

Let  u0(x,*)  denote  the  solution  of  (7.1),  with  <  replac¬ 
ing  V  q ,  and  let  u1(l),  u2 ( A  ) ,  ....  un_.(\),  v^ ( v  )  be  the 
corresponding  sequence  obtained  from  the  finite  difference 
solution.  Let  c(A)  denote  the  number  of  sign  changes  in  this 
sequence.  Choose  A,  A'  close  to  =  'k*  such  that 

\  <  \  q  <  \  .  By  Case  1 ,  c  ( A  )  =  k- 1  and  c  (  v. '  )  =  k .  Since 
Uq(1,Ao)  *  0,  the  number  of  zeros  in  (0,1)  is  the  same  for 
uQ(x,v),  Uq(x,Aq),  Uq(x,A').  Also,  this  number  of  zeros  equals 
the  number  of  sign  changes  in  the  sequences 

u1(A),  u2(A),  ...,  un_1(A)  and  UjfA'),  u2  ( A  )  ,  ...,  un_1(A'). 

Since  c ( A ' )  =  c(A)  +  1,  the  extra  sign  change  must  come  from  a 

has  k- 1  zeros 


/.v. 


change  in  sign  of  vn  .  This  shows  that  u Q  \  ^  ) 


in  (0,1),  which  proves  the  theorem  in  Case  3. 

Case  4 .  uQ(l)  =0,  v  =  0. 

Again  is  an  eigenvalue,  say  Ao  =  'k‘  In  case'  we 

have  a  Dirichlet  boundary  condition:  a1  u{l)  +  1 1  u'  ( 1 )  -  u(l)  ; 
v  =  Uq(1),  and  a  =  0  by  definition.  We  must  show  that 
Uq(x,\q)  has  k-1  zeros  in  (0,1). 

As  in  the  previous  case,  choose  A,  A'  close  to  \q,  with 


A 

A 

>• 

o 

A 

7* 

Then 

V  ( A  ) 

=  Uq  ( 1 ,  A. )  *  0  and  v  (  A  '  )  = 

u0 

(1. 

V'  ) 

x  0 

By 

Case  1 , 

uQ ( x , A ) 

has 

k-1  zeros  in  (0,1)  and 

u0 

(X, 

' ' ) 

ha 

k 

zeros  in 

(0,1) . 

The 

extra  zero  must  arise  from 

u0 

(x, 

Ko') 

,  a 

x  =  1.  Therefore  uQ(x,A0)  has  k-1  zeros  in  (0,1).  This 
concludes  the  proof  of  the  theorem. 


8.  Conclusions. 


A  closure  exists  for  any  numerical  process  dealing  with  a 
matrix  problem  which  comes  from  discretizing  an  ordinary  aifferen 
tial  equation.  We  have  shown  that  the  closure  of  the  Sturm 
sequence  algorithm  for  finding  eigenvalues  leads  to  a  certain 
shooting  method.  Other  numerical  processes,  such  as  vector  iter¬ 
ation  methods,  transformation  methods  and  polynomial  iteration 
techniques  also  have  closures.  We  have  mentioned  some  closures 
(such  as  the  double  sweep  method)  of  algorithms  for  solving 
systems  of  linear  equations.  A  closure  of  an  algorithm  is  impor¬ 
tant,  because  it  leads  to  a  new  method  for  solving  the  continuous 
problem.  The  new  method  is  often  superior  to  the  finite  differ¬ 
ence  or  finite  element  method,  because  it  has  greater  flexibility 
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